package scu.maqiang.fes;

import java.util.Arrays;

import scu.maqiang.mesh.MatrixFunc;
import scu.maqiang.mesh.Mesh2T6;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.CG;
import scu.maqiang.numeric.Constants;
import scu.maqiang.numeric.DoubleIntegerPair;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.ParamCheck;
import scu.maqiang.numeric.SRMatrix;

public class FES2T61 extends FES {
	public FES2T61(Mesh2T6 mesh) {
		super(mesh, 1);
		eleStiff = FES2T61::elementStiff;
		eleMass = FES2T61::elementMass;
		eleSource = FES2T61::elementSource;
	}

	public double getValue(double[] u, double[] xy, int[] eleIdx, double[] eleVal) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
	    double[] uu = new double[6];
	    double value = 0;
	    double[][] coord = new double[6][2];
	    for(int i = 0, n = mesh.getNt(); i < n; i++) {
	    	int[] ele = mesh.getElement(i);
	    	mesh.getCoordPerElement(ele, coord);
	    	CG.T6ShapeFunction(coord, xy, eleVal);
	    	//System.out.println(Arrays.toString(eleVal));
	    	if (Math.abs(MVO.sum(eleVal) - 1.0) < 1.0e-10) {
	    		for(int j = 0; j < 6; j++) {
	    			uu[j] = u[ele[j]];
	    			eleIdx[j]= ele[j];
	    			value += uu[j] * eleVal[j];
	    		}
	    		//System.out.println(value);
	    		return value;
	    	}
	    }
	    throw new IllegalArgumentException("The Point ( " + Arrays.toString(xy) + " ) is not in the domain!");
	}
	
	public void getValue(double[] u, double[] xy, DoubleIntegerPair data, int[] eleIdx, double[] eleVal) {
		ParamCheck.checkEqual(u.length, this.ndof);
		ParamCheck.checkEqual(eleIdx.length, this.nPerEle);
		ParamCheck.checkEqual(eleVal.length, this.nPerEle);
	    double value = 0;
		int[] ele = mesh.getElement(data.second);
		mesh.getCoordPerElement(ele, coordRegion);
		CG.T6ShapeFunction(coordRegion, xy, eleVal);
		if (Math.abs(MVO.sum(eleVal) - 1.0) < 1.0e-10) {
			data.first = MVO.dot_product(u, ele, eleVal);
			return;
		} else {
		    for(int i = 0, n = mesh.getNt(); i < n; i++) {
		    	ele = mesh.getElement(i);
		    	mesh.getCoordPerElement(ele, coordRegion);
		    	CG.T6ShapeFunction(coordRegion, xy, eleVal);
		    	//System.out.println(Arrays.toString(eleVal));
		    	if (Math.abs(MVO.sum(eleVal) - 1.0) < 1.0e-10) {
		    		for(int j = 0; j < 6; j++) {
		    			eleIdx[j]= ele[j];
		    			value += u[ele[j]] * eleVal[j];
		    		}
		    		data.first = value;
		    		data.second = i;
		    		return;
		    	}
		    }
		    throw new IllegalArgumentException("The Point ( " + Arrays.toString(xy) + " ) is not in the domain!");
		}

	   
	}

	/**
	 * ���㵥Ԫ�նȾ���
	 * 
	 * @param xy
	 * @param kxx
	 * @param tp
	 * @param eleMatrix
	 */
	public static void elementStiff(int[] ele, double[][] xy, double[][] kxx, BVPType tp, double[][] eleMatrix) {

		double[] x = MVO.col(xy, 0);
		double[] N = new double[6];
		double[][] gradN = new double[6][2];
		double[] L = new double[3];
		double area, weight;
		MVO.fill(eleMatrix, 0.0);
		switch (tp) {
		case COMMON:
			for (int i = 0; i < 4; i++) {
				L = Constants.TP3[i].clone();
				area = CG.T6ShapeFunction(xy, L, N, gradN);
				weight = area * MVO.dot_product(kxx[0], N) * Constants.TW3[i];
				for (int si = 0; si < 6; si++) {
					for (int sj = 0; sj < 6; sj++) {
						eleMatrix[si][sj] += (gradN[si][0] * gradN[sj][0] + gradN[si][1] * gradN[sj][1]) * weight;
					}
				}
			}
			break;
		case AXISYMMETRIC:
			for (int i = 0; i < 4; i++) {
				L = Constants.TP3[i].clone();
				area = CG.T6ShapeFunction(xy, L, N, gradN);
				weight = area * MVO.dot_product(kxx[0], N) * MVO.dot_product(x, N) * Constants.TW3[i];
				for (int si = 0; si < 6; si++) {
					for (int sj = 0; sj < 6; sj++) {
						eleMatrix[si][sj] += (gradN[si][0] * gradN[sj][0] + gradN[si][1] * gradN[sj][1]) * weight;
					}
				}
			}
			break;
		case CONSTITUTE_MATRIX_COMMON:
			for (int m = 0; m < 3; ++m) {
				L = Constants.TP2[m].clone();
				area = CG.T6ShapeFunction(xy, L, N, gradN);
				double[][] temp = MVO.matmul(gradN, kxx);
				weight = area * Constants.TW2[m];
				for (int i = 0; i < 6; i++) {
					for (int j = 0; j < 6; j++) {
						for (int k = 0; k < 2; k++) {
							eleMatrix[i][j] += temp[i][k] * gradN[j][k] * weight;
						}
					}
				}
			}
			break;
		case CONSTITUTE_MATRIX_AXIS:
			for (int m = 0; m < 3; m++) {
				L = Constants.TP2[m].clone();
				area = CG.T6ShapeFunction(xy, L, N, gradN);
				double[][] temp = MVO.matmul(gradN, kxx);
				weight = area * Constants.TW2[m] * MVO.dot_product(x, N);
				for (int i = 0; i < 6; i++) {
					for (int j = 0; j < 6; j++) {
						for (int k = 0; k < 2; k++) {
							eleMatrix[i][j] += temp[i][k] * gradN[j][k] * weight;
						}
					}
				}
			}
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
		
		//System.out.println(MVO.toString(eleMatrix));
		//System.exit(0);
	}

	
	public static void elementMass(int[] ele, double[][] xy, double[][] coef, BVPType tp, double[][] eleMatrix) {
		double[] N = new double[6];
		double[][] gradN = new double[6][2];
		MVO.fill(eleMatrix, 0.0);

		for (int i = 0; i < 3; i++) {
			double[] L = Constants.TP2[i].clone();
			double area = CG.T6ShapeFunction(xy, L, N, gradN);
			double weight = area * MVO.dot_product(coef[0], N) * Constants.TW2[i];
			for (int si = 0; si < 6; si++) {
				for (int sj = 0; sj < 6; sj++) {
					eleMatrix[si][sj] += N[si] * N[sj] * weight;
				}
			}
		}
		if (tp == BVPType.CONSIST_MASS) {
			return;
		} else if (tp == BVPType.LUMP_MASS) {
			double[] sumMatrix = MVO.sumRow(eleMatrix);
			MVO.fill(eleMatrix, 0.0);
			for (int i = 0; i < 6; i++) {
				eleMatrix[i][i] = sumMatrix[i];
			}
		}
	}


	public static void elementSource(int[] ele, double[][] xy, double[][] f, BVPType tp, double[] eleF) {
		double[] x = MVO.col(xy, 0);
		double[] N = new double[6];
		double[][] gradN = new double[6][2];
		double[] L = new double[3];
		double area, weight;
		Arrays.fill(eleF, 0.0);
		switch (tp) {
		case COMMON:
		case CONSTITUTE_MATRIX_COMMON:
			for (int i = 0; i < 3; i++) {
				L = Constants.TP2[i].clone();
				area = CG.T6ShapeFunction(xy, L, N, gradN);
				weight = area * MVO.dot_product(f[0], N) * Constants.TW2[i];
				for (int si = 0; si < 6; si++) {
					eleF[si] += N[si] * weight;
				}
			}
			break;
		case AXISYMMETRIC:
		case CONSTITUTE_MATRIX_AXIS:
			for (int i = 0; i < 3; i++) {
				L = Constants.TP2[i].clone();
				area = CG.T6ShapeFunction(xy, L, N, gradN);
				weight = area * MVO.dot_product(f[0], N) * MVO.dot_product(x, N) * Constants.TW2[i];
				for (int si = 0; si < 6; si++) {
					eleF[si] += N[si] * weight;
				}
			}
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
	}
	
	public void computeGradient(double[] u, double[][] gradU) {
		int len = u.length;
		ParamCheck.checkEqual(len, this.ndof);
		ParamCheck.checkEqual(len, gradU[0].length);
		ParamCheck.checkEqual(len, gradU[1].length);
		
		double[] angle = new double[6];
		double[][] coord = new double[6][2];
		double[][] grad = new double[3][2];
		double[][] eleGradU = new double[6][2];
		double[] angleCount = new double[this.GetNdof()];
		double[] eleU = new double[6];
		for(int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coord);
			for(int j = 0; j < 6; j++) {
				eleU[j] = u[ele[j]];
			}
			CG.gradient2DT3(coord, grad);
			CG.gradientValue2DT6(grad, eleU, eleGradU);
			CG.angle2DT6(coord, angle);
			for(int j = 0; j < 6; j++) {
				gradU[0][ele[j]] += eleGradU[j][0] * angle[j];
				gradU[1][ele[j]] += eleGradU[j][1] * angle[j];
				angleCount[ele[j]] += angle[j];
			}
			//System.out.println(Arrays.toString(MVO.col(eleGradU, 1)));
		}
		
		//System.out.println(MVO.toString(gradU[0]));
		for(int i = 0, n = mesh.getNv(); i < n; i++) {
			gradU[0][i] /= angleCount[i];
			gradU[1][i] /= angleCount[i];
		}
	}
	
	public void computeSecondDerivative(double[] u, double[][] sdU) {
		int len = u.length;
		ParamCheck.checkEqual(len, this.ndof);

		double[] angle = new double[6];
		double[][] coord = new double[6][2];
		double[][] grad = new double[3][2];
		double[] eleSDU = new double[3];
		double[] angleCount = new double[this.GetNdof()];
		double[] eleU = new double[6];
		for(int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coord);
			for(int j = 0; j < 6; j++) {
				eleU[j] = u[ele[j]];
			}
			CG.gradient2DT3(coord, grad);
			CG.secondDerivative2DT6(grad, eleU, eleSDU);
			CG.angle2DT6(coord, angle);
			for(int j = 0; j < 6; j++) {
				sdU[0][ele[j]] += eleSDU[0] * angle[j];
				sdU[1][ele[j]] += eleSDU[1] * angle[j];
				//sdU[2][ele.idx[j]] += eleSDU[2] * angle[j];
				sdU[2][ele[j]] += eleSDU[2] * angle[j];
				angleCount[ele[j]] += angle[j];
			}
		}
		for(int i = 0, n = mesh.getNv(); i < n; i++) {
			sdU[0][i] /= angleCount[i];
			sdU[1][i] /= angleCount[i];
			sdU[2][i] /= angleCount[i];
		}
		
//		for(int i = 0; i < mesh.getNt(); i++) {
//			ele = mesh.getElement(i);
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdU[0][ele.idx[j]];
//			}
//			sdU[0][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdU[0][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdU[0][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//			
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdU[1][ele.idx[j]];
//			}
//			sdU[1][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdU[1][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdU[1][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//			
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdU[2][ele.idx[j]];
//			}
//			sdU[2][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdU[2][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdU[2][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//		}
	}
	
	public void computeFirstAndSecondDerivative(double[] u, double[][] fdU, double[][] sdU) {
		int len = u.length;
		ParamCheck.checkEqual(len, this.ndof);

		double[] angle = new double[6];
		double[][] coord = new double[6][2];
		double[][] grad = new double[3][2];
		double[][] eleGradU = new double[6][2];
		double[] eleSDU = new double[3];
		double[] angleCount = new double[this.GetNdof()];
		double[] eleU = new double[6];
		for(int i = 0, n = mesh.getNt(); i < n; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coord);
			for(int j = 0; j < 6; j++) {
				eleU[j] = u[ele[j]];
			}
			CG.gradient2DT3(coord, grad);
			CG.gradientValue2DT6(grad, eleU, eleGradU);
			CG.secondDerivative2DT6(grad, eleU, eleSDU);
			CG.angle2DT6(coord, angle);
			for(int j = 0; j < 6; j++) {
				fdU[0][ele[j]] += eleGradU[j][0] * angle[j];
				fdU[1][ele[j]] += eleGradU[j][1] * angle[j];
				sdU[0][ele[j]] += eleSDU[0] * angle[j];
				sdU[1][ele[j]] += eleSDU[1] * angle[j];
				sdU[2][ele[j]] += eleSDU[2] * angle[j];
				angleCount[ele[j]] += angle[j];
			}
		}
		for(int i = 0, n = mesh.getNv(); i < n; i++) {
			fdU[0][i] /= angleCount[i];
			fdU[1][i] /= angleCount[i];
			sdU[0][i] /= angleCount[i];
			sdU[1][i] /= angleCount[i];
			sdU[2][i] /= angleCount[i];
		}
		
		for(int i = 0; i < mesh.getNt(); i++) {
			int[] ele = mesh.getElement(i);
			for(int j = 0; j < 6; j++) {
				eleU[j] = sdU[0][ele[j]];
			}
			sdU[0][ele[3]] = 0.5 * (eleU[0] + eleU[1]);
			sdU[0][ele[4]] = 0.5 * (eleU[1] + eleU[2]);
			sdU[0][ele[5]] = 0.5 * (eleU[0] + eleU[2]);
			
			for(int j = 0; j < 6; j++) {
				eleU[j] = sdU[1][ele[j]];
			}
			sdU[1][ele[3]] = 0.5 * (eleU[0] + eleU[1]);
			sdU[1][ele[4]] = 0.5 * (eleU[1] + eleU[2]);
			sdU[1][ele[5]] = 0.5 * (eleU[0] + eleU[2]);
			
			for(int j = 0; j < 6; j++) {
				eleU[j] = sdU[2][ele[j]];
			}
			sdU[2][ele[3]] = 0.5 * (eleU[0] + eleU[1]);
			sdU[2][ele[4]] = 0.5 * (eleU[1] + eleU[2]);
			sdU[2][ele[5]] = 0.5 * (eleU[0] + eleU[2]);
		}
	}
	
	
//subroutine computeSecondDerivativeAngleT61(this, T, dxxT, dxyT, dyxT, dyyT)
//    implicit none
//    class(FES2T61), intent(in) :: this
//    real*8, intent(in) :: T(:)
//    real*8, intent(out) :: dxxT(:), dxyT(:), dyxT(:), dyyT(:)
//
//    real*8 :: x(6), y(6), dx(6), dy(6), angle(6)
//    real*8 :: uu(6), dxxU, dxyU, dyxU, dyyU
//    integer :: element(6)
//    type(Node2D) :: node(6)
//    type(Mesh2T6), pointer :: mesh2D
//    real*8 :: area, dxValue, dyValue
//    real*8, allocatable :: angleCount(:)
//    integer :: i
//    mesh2D => this.mesh2D
//    
//    allocate(angleCount(this%nDof))
//    angleCount = 0.0D0
//    dxxU = 0.0D0
//    dxyU = 0.0D0
//    dyxU = 0.0D0
//    dyyU = 0.0D0
//    dxxT = 0.0D0
//    dxyT = 0.0D0
//    dyxT = 0.0D0
//    dyyT = 0.0D0
//    do i = 1, this%nElement
//        element = mesh2D%elements(i)%idx
//        node = mesh2D%nodes(element)
//        x = node%xy(1)
//        y = node%xy(2)
//        uu = T(element)
//        call gradient2DP1(x, y, dx, dy, area)
//        call secondDerivativeValueP2(x, y, dx, dy, uu, dxxU, dxyU, dyxU, dyyU)
//        !call gradientValueP2(x, y, dx, dy, uu, dxU, dyU)
//        angle = angleTriangleP2(x, y)
//        !!write(*, *) "angle", angle, sum(angle)
//
//        dxxT(element) = dxxT(element) + dxxU * angle
//        dxyT(element) = dxyT(element) + dxyU * angle
//        dyxT(element) = dyxT(element) + dyxU * angle
//        dyyT(element) = dyyT(element) + dyyU * angle
//        angleCount(element) = angleCount(element) + angle
//    end do
//    dxxT = dxxT / angleCount
//    dxyT = dxyT / angleCount
//    dyxT = dyxT / angleCount
//    dyyT = dyyT / angleCount
//
//
//    !!write(*, *)"angleCount", angleCount
//    deallocate(angleCount)
//
//
//    do i = 1, this%nElement
//        element = mesh2D%elements(i)%idx
//        node = mesh2D%nodes(element)
//        x = node%xy(1)
//        y = node%xy(2)
//        uu = dxxT(element)
//        uu(4) = 0.5 * (uu(1) + uu(2))
//        uu(5) = 0.5 * (uu(2) + uu(3))
//        uu(6) = 0.5 * (uu(1) + uu(3))
//        dxxT(element) = uu
//
//        uu = dyyT(element)
//        uu(4) = 0.5 * (uu(1) + uu(2))
//        uu(5) = 0.5 * (uu(2) + uu(3))
//        uu(6) = 0.5 * (uu(1) + uu(3))
//        dyyT(element) = uu
//
//        uu = dxyT(element)
//        uu(4) = 0.5 * (uu(1) + uu(2))
//        uu(5) = 0.5 * (uu(2) + uu(3))
//        uu(6) = 0.5 * (uu(1) + uu(3))
//        dxyT(element) = uu
//    end do
//    dyxT = dxyT
//end subroutine
	
	
//	public void computeFirstAndSecondDerivative(double[] u, double[][] fdU, double[][] sdU) {
//		int len = u.length;
//		ParamCheck.checkEqual(len, this.ndof);
//
//		double[] angle = new double[6];
//		double[][] coord = new double[6][2];
//		double[][] grad = new double[3][2];
//		double[][] eleGradU = new double[6][2];
//		double[] eleSDU = new double[4];
//		double[] angleCount = new double[this.GetNdof()];
//		Element ele = null;
//		double[] eleU = new double[6];
//		for(int i = 0, n = mesh.getNt(); i < n; i++) {
//			ele = mesh.getElement(i);
//			mesh.getCoordPerElement(ele, coord);
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = u[ele.idx[j]];
//			}
//			CG.gradient2DT3(coord, grad);
//			CG.gradientValue2DT6(grad, eleU, eleGradU);
//			CG.secondDerivative2DT6(grad, eleU, eleSDU);
//			CG.angle2DT6(coord, angle);
//			for(int j = 0; j < 6; j++) {
//				fdU[0][ele.idx[j]] += eleGradU[j][0] * angle[j];
//				fdU[1][ele.idx[j]] += eleGradU[j][1] * angle[j];
//				sdU[0][ele.idx[j]] += eleSDU[0] * angle[j];
//				sdU[1][ele.idx[j]] += eleSDU[1] * angle[j];
//				sdU[2][ele.idx[j]] += eleSDU[2] * angle[j];
//				angleCount[ele.idx[j]] += angle[j];
//			}
//		}
//		for(int i = 0, n = mesh.getNv(); i < n; i++) {
//			fdU[0][i] /= angleCount[i];
//			fdU[1][i] /= angleCount[i];
//			sdU[0][i] /= angleCount[i];
//			sdU[1][i] /= angleCount[i];
//			sdU[2][i] /= angleCount[i];
//		}
//		
//		for(int i = 0; i < mesh.getNt(); i++) {
//			ele = mesh.getElement(i);
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdU[0][ele.idx[j]];
//			}
//			sdU[0][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdU[0][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdU[0][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//			
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdU[1][ele.idx[j]];
//			}
//			sdU[1][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdU[1][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdU[1][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//			
//			for(int j = 0; j < 6; j++) {
//				eleU[j] = sdU[2][ele.idx[j]];
//			}
//			sdU[2][ele.idx[3]] = 0.5 * (eleU[0] + eleU[1]);
//			sdU[2][ele.idx[4]] = 0.5 * (eleU[1] + eleU[2]);
//			sdU[2][ele.idx[5]] = 0.5 * (eleU[0] + eleU[2]);
//		}
//	}
}
